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The phase diagram of the massive chiral Gross-Neveu model (the massive Nambu-Jona-Lasinio 
model in 1+1 dimensions) is constructed. In the large A'' limit, the Hartree-Fock approach can 
be used. We find numerically a chiral crystal phase separated from a massive Fermi gas phase by 
a 1st order transition. Using perturbation theory, we also construct the critical sheet where the 
homogeneous phase becomes unstable in a 2nd order transition. A tricritical curve is located. The 
phase diagram is mapped out as a function of fermion mass, chemical potential and temperature 
and compared with the one of the discrete chiral Gross-Neveu model. As a by-product, we illustrate 
the crystal structure of matter at zero temperature for various densities and fermion masses. 
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I. INTRODUCTION 

To map out the phase diagram of hot and dense mat- 
ter has been a major goal of strong interaction physics 
during the last decades, both experimentally and theo- 
retically. As is often the case, these efforts have been 
accompanied by studies of drastically simplified, solvable 
model problems to sharpen the theoretical tools and get 
guidance for more realistic cases. Among the few known 
field theories which are both solvable and possess a non- 
trivial phase structure, fermionic large N models in 1-1-1 
dimensions like the 't Hooft model [l| or Gross-Neveu 
models [2] are perhaps most instructive, as they share 
a number of properties with quantum chromodynamics 
(for a pedagogical review, see Ref. [3]). Given that these 
models have been formulated back in 1974 already, it 
is surprising that their phase diagrams as a function of 
temperature, chemical potential and fermion mass have 
not yet been fully established. As far as we can tell, 
the reason is not that these phase diagrams were consid- 
ered to be uninteresting. Rather, this situation reflects 
a shortcoming of the first round of theoretical investiga- 
tions during the 80's and 90's with methods too crude to 
expose the full, rich phase structure. As a consequence, 
there has been renewed interest recently in this topic with 
results which also have some bearing on low dimensional 
condensed matter systems, and the original phase dia- 
grams are still in the process of revision right now. For 
an update on the current state of the art, see Refs. [J, Q 
and references therein. 

In the present paper, we focus on the phase structure 
of the massive, chiral Gross-Neveu (GN) model at finite 
temperature and chemical potential. This model is noth- 
ing but the 1-1-1 dimensional Nambu-Jona-Lasinio model 
[6| with N fermion flavors and a bare mass term explicitly 
breaking the U(l)(g)U(l) chiral symmetry. Its Lagrangian 



reads 



C = V5(i^ - mo)i> + ^ [(^V)' + (^iTsV')'] (1) 
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where flavor indices are suppressed as usual, i.e., 'ip4' — 
J2k=i '^k'4'k etc. We are only interested in the 't Hooft 
limit {N — > OO, Ng'^ = const.) in which classic no-go 
theorems can be bypassed and breakdown of continuous 
symmetries becomes possible in 1-1-1 dimensions. In spite 
of the fact that semiclassical methods for solvingsuch 
models have been developed in the 70's already [2, HIj 
the full phase diagram of the simple field theory with 
Lagrangian ([1]) is still largely unknown. Consider first 
the chiral limit (too = 0) of the model. If one constrains 
the condensates {4>ip) , {ip^'j^^') to be spatially constant, 
the resulting phase diagram is identical to the one from 
the simpler GN model variant with discrete chiral sym- 
metry and scalar-scalar coupling (V'V')^ only. One finds 
two phases in the (/x, T) plane, a massless and a massive 
Fermi gas, separated by 1st and 2nd order transitions [g] . 
As soon as one allows for spatially inhomogeneous con- 
densates, the system takes advantage of the Peierls effect 
9] and opens a gap at the Fermi surface. This results in 
a solitonic chiral crystal phase. The first crystalline solu- 
tion of the Hartree-Fock (HF) problem which was found 
is the "chiral spiral" with helical order parameter and a 
strikingly different phase diagram [J, |T3| . As pointed out 
in Ref. [llj , these results can also be understood readily 
in terms of bosonization. Very recently however, they 
have been challenged by a more sophisticated candidate 
for the complex order parameter in the form of a chi- 
rally twisted crystal, using powerful resolvent methods 
to generate self-consistent solutions in closed analytical 
form [^, [ij] . The implications for the phase diagram have 
not yet been fully worked out but promise an even richer 
structure of the solitonic crystal phase than previously 
thought. 

Turning to the massive chiral GN model (mo > 0), we 
first should like to remind the reader that the bare pa- 
rameters g^,mo in Eq. ([T]) together with the UV cutoff 
A/2 get replaced by two physical, renormalization group 
invariant parameters m and 7 in the process of regular- 



ization and renornialization [J,|l4|- Here, m is the phys- 
ical fermion mass in the vacuum (set equal to 1 without 
loss of generality throughout this paper) and 7 the "con- 
finement parameter" measuring the explicit violation of 
chiral symmetry, 



7Vg= 



= 7- 



in A, 

m 



7 := 



TT ?71o 



(2) 



The following bits and pieces are known about the phase 
structure of the massive model. The phase diagram as- 
suming x-independent condensates only [ij| is once again 
indistinguishable from that of the massive discrete chiral 
GN model, but this is an artefact of the assumption of 
homogeneity Q. As far as inhomogeneous condensates 
are concerned, it is useful to start from the low density, 
low temperature limit governed by the isolated baryons of 
the model. Baryons of the massive chiral GN model were 
first studied near the chiral limit by means of variational 
techniques [l5| and subsequently via the derivative ex- 
pansion [1^. They turn out to be closely related to the 
sine-Gordon kink. A recent numerical HF calculation, 
supplemented by analytical asymptotic expansions, has 
been able to follow the baryon mass and structure to arbi- 
trary 7 [17.] . This was actually done in preparation of the 
present study. Unlike in the discrete chiral GN model, 
the self-consistent baryon potentials found were not re- 
flectionless, a serious obstacle for a full analytical solu- 
tion. Aside from individual baryons relevant to the base 
line at T = of the (7, /i, T) phase diagram, the vicinity 
of the tricritical point (7 = 0, u = 0,T = c^/tt) has also 
been explored in some detail [l8|. The phase structure 
was deduced from a microscopic Ginzburg-Landau (GL) 
approach, based once again on the derivative expansion. 
In this work, both first and second order critical lines 
between homogeneous and inhomogeneous phases were 
identified. As a result, one already starts to see that the 
GN models with (broken) discrete and continuous chi- 
ral symmetry have totally different phase diagrams, as is 
indeed expected on the basis of universality arguments. 

In the present paper, we report on a solution of the HF 
problem at finite T, /i for a whole range of 7 values and 
construct a first candidate for the full phase diagram of 
the massive chiral GN model. It is not yet known what 
impact the more general, chirally twisted soliton crystals 
of the massless model discovered in Ref. [13] would have 
on the massive model, and we cannot contribute anything 
to this question. Our aim here is to extend the calcula- 
tions of Ref. 18] near the tricritical point to a significant 
portion of (/i, T, 7) space, so that a 3d plot of the phase 
diagram can be drawn and compared with the one from 
the discrete chiral GN model. We think that such an 
undertaking is worthwhile in the present situation, but 
should be followed up by efforts to identify alternative 
chiral crystal structures which might be thermodynam- 
ically more stable [13|, or by further attempts to arrive 
at a full analytical solution as in the case of the discrete 
chiral GN model [1] . 

The remaining paper is organized as follows. Sec. Ullis 



devoted to the HF calculation at zero temperature. We 
explain the general numerical procedure (jll A[) . discuss 
analytically the low and high density asymptotics (jIIB[) 
and present selected numerical results (jIICp . Sec. IIIII 
contains all the material about finite temperature and 
the phase diagram. We briefly outline the thermal HF 
approach to the grand canonical potential (jIII Ap and re- 
call previous results from GL theory piIBp . In Sec. lIIICl 
we describe in detail how we have obtained the pertur- 
bative, 2nd order critical sheet. The non-perturbative 
1st order sheet represents the most difficult part of our 
analysis, since we can only determine it numerically at 
present. This is presented in Sec. IIIIDJ along with the fi- 
nal results. In the concluding Sec. llVi we summarize our 
findings, compare the phase diagram with other related 
phase diagrams and identify areas where more work is 
needed. 



II. HARTREE-FOCK CALCULATION OF 
DENSE MATTER AT T = 

A. Setup of the numerical calculations 

The HF calculation in the chiral GN model starts from 
the Dirac Hamiltonian 

1 d 



H 



75- 



dx 



fS{x)+i-i^P{x) 



(3) 



with scalar and pseudoscalar potentials 5, P to be de- 
termined self-consistently. In Ref. [i3], a numerical HF 
study including the Dirac sea has been used to construct 
the baryons of this model. For technical reasons, the cal- 
culation was done in a finite interval of length L with 
antiperiodic boundary conditions for the fermions, using 
a basis of free, massive spinors in discretized momentum 
space. Now assume that 5*, P are periodic with spatial 
period a. We can actually reduce the HF calculation 
for such a crystal to the one for a single baryon per- 
formed in [13 . We enclose the crystal in a box of length 
L = Na containing N periods and impose again antiperi- 
odic boundary conditions on the fermion single particle 
wave functions in this large interval. 



V'(L) = -V'(O). 



(4) 



According to the Bloch theorem, the eigenspinors of H 
are of the form 



-0(a;) = (j){x)e^' 



[X + a) = mix) 



(5) 



The boundary condition (JH) discretizes the Bloch mo- 
menta. 
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P„ = - n+- 



1 



{n e Z). 



(6) 



For a single period, e.g., the interval [0, a], this implies 
quasi-periodic boundary conditions. 



^(a) = e''^-V(0), 



(7) 



where the N discrete values of /3y parametrize the 7V-th 
roots of (—1), 



Pu 



2tt 



0,1, 



,N-l. 



(8) 



Hence, to get the spectrum of H with a periodic po- 
tential, all we have to do is compute the spectrum for 
a single "baryon" in an interval of length a with quasi- 
periodic boundary conditions along the lines of Ref. [l3] , 
repeat the calculation N times (for all possible values 
of the phase P^) and collect the spectra. This enables 
us to take over the calculational method literally from 
Ref. 113. We also stick to the conditions 



S{x) = S{~x), P{x) = "P(-x), 



(9) 



reflecting the difference between scalar and pseudoscalar 
potentials if parity is unbroken. To evaluate the energy 
density of the crystal at T = 0, we once again combine 
a numerical diagonalization with perturbation theory for 
states deep down in the Dirac sea. The technical details 
like vacuum subtraction, double counting correction and 
renormalization are identical to those given in Ref. |17| 
and need not be repeated here. 

A key element of the HF approach is self-consistency 
of the potentials S, P. As explained in Ref. [l3| this 
can be achieved by minimizing the HF energy at fixed 
fermion number with respect to the potentials, provided 
one varies the potentials without any bias. In the present 
work, we assume periodicity, expand S and P into Fourier 
series. 



S{x) = J2 ^^e'^-^^/-^, P{x) = i^ P^e'2- 



£x/a 



(10) 



and minimize the HF energy with respect to the Fourier 
coefficients Se, Pi and the spatial period a, using a stan- 
dard conjugate gradient algorithm. The only other bias 
put in aside from periodicity are the symmetry relations 
^. If the true self-consistent potential A = S—iP would 
not be strictly periodic but carry a chiral twist, 



A{x + a)=e'^''<'A{x), 



(11) 



as proposed in a recent study of the massless chiral GN 
model [H, [l3| , our calculation might still be useful as a 
variational calculation, but we could miss the true self- 
consistent potential. Note however that there is so far 
no claim of non-periodic potentials in the massive model 
considered in the present work. 



B. Low and high density limits 

In the limits of low and high fermion density, the 
ground state energy can be calculated analytically. If 
the valence band is completely filled (as is indeed found 



in the full HF calculation) , the spatially averaged baryon 
density per flavor is related to the period a via 



1 



P = 






(12) 



The last equation defines the Fermi momentum pf. At 
very low density, we expect the energy density to be de- 
termined by the baryon mass. 



Shf — Sv 



MbP 



(13) 



The baryon mass is known already from Ref. [17| . At 
high density on the other hand, we can use perturbation 
theory to predict the asymptotic behaviour of the en- 
ergy density. This is a simple generalization of a similar 
calculation done in Ref. [l9| for the massive GN model 
with broken discrete chiral symmetry, cf. Eqs. (67)-(74) 
of that paper. It is sufficient to keep the Fourier am- 
plitudes So, Si and Pi for this purpose. Standard 2nd 
order perturbation theory then yields the single particle 
energies 

/ ^2 , iSi+fPir , (Si-fPi)' 

E,,,p = 77Sgn(p) [p+TT 



2p 2{p + pf) 2(p-p/) 



where 77 = ±1 and 



/ ^ ^^ri^sgn(p) • 



(14) 



(15) 



Along the lines of Ref. [l^ , we find for the perturbative 
ground state energy 



AV p/ , ^o^ 



75*0 



£hf = -^ + ^ + :r^[7 + ln(2p/)] - -^ (16) 

on Zn ZTT TT 

+ iL[2^_i + in(y2)] + [7 + ln(4p;)] 

47r TT 

where we have set 

Si=X + y/2, Pi=X-y/2. (17) 

Minimizing fuF with respect to Sq, X and y yields 

7 



X = 0, Sn = 



(18) 



j + H2pf) 

and the equation 

y[27-fln(y2)]=0 (19) 

with the solutions y = (homogeneous condensate) and 

y = ±e-\ (20) 

The self-consistent potential A — S* — iP for the non- 
trivial solution y — e~'^ is inhomogeneous. 



A(x) 



7 



7 -I- \n{2pf) 



exp{— 2ip/a: ~ 7} . (21) 



(The other sign of y merely corresponds to a translation 
of the crystal by a/2.) The ground state energy (fTB|) at 
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FIG. 1: Self-consistent scalar HF potential S{x) at T = 
0,Pf = 0.2 and 7 = 0.2,0.6, 1.0 (from bottom to top), showing 
well resolved baryons. Here and in Figs. [2H4] only one spatial 
period of the periodic potentials is shown. 
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FIG. 2: One period of self-consistent pseudoscalar HF po- 
tential P{x) at T = 0,pf = 0.2 and 7 = 0.2,0.6,1.0 (with 
decreasing amplitude). 



the minimuni is indeed lower than the one of the homo- 
geneous solution, 



fHF(y^O)-£HF(y = 0) 



1 

6 

An 



-27 



(22) 



At 7 = this agrees with the result for the chiral spiral 
[lO|. Finally we write down the ground state energy for 
large pf, relative to the vacuum. It has the asymptotic 
behavior 



p} 

^HF — f vac ~ 7;— 



7 



27r 27r(7-|-ln2p/) 



+ ^(1 



27- 



.-27\ 



(23) 



Eqs. (I13|) for Pf -^ and ([23]) for pf -^ cx3 are the main 
results of this section, ready to be compared to full nu- 
merical results below. 



C. Numerical results 

We vary with respect to the Fourier components Si, P( 
which, owing to Eqs. ([9]) and (fTO| . are real and satisfy 
S^i — Si, P-i = —Pi. The actual calculations were done 
as follows. We use N = L/a = 8, i.e., perform single 
baryon computations with 8 different boundary condi- 
tions. In the sum over single particle energies we now 
have to subtract numbers of 0(100 000), as compared 
to 0(10 000) for a single baryon. To keep computations 
feasible with MAPLE, we had to compromise on the size 
of the momentum space basis and choose the smallest 
size which gave sufficient precision in the tests, TV = 50 
(corresponding to 201x201 matrices). The total number 
of single particle states computed by diagonalization is 
therefore 8 x 201 = 1608. We kept all Fourier modes 



of S, P up to ^max — 6, SO that 14 real parameters had 
to be varied in total (the period a, 5*0, and {Si,Pe} for 
£ — 1...6). To test our MAPLE code, we computed the 
energy density of crystals where the analytic solution is 
known (chiral spiral, massless and massive GN models). 
In all of these cases the energy density was reproduced 
correctly to 7 significant digits. The other uncertainty 
comes from the minimization procedure. It was found 
that after only 20 conjugate gradient steps, the results as 
shown in the figures below did not change anymore sig- 
nificantly. Under these conditions, all calculations could 
still be done using MAPLE on high-end PC's, without 
need to switch to compiled programming languages. 

We now turn to the results of the T = computations. 
Just as in the discrete chiral GN model, we found that 
it is always energetically advantageous to let the Fermi 
surface coincide with the lower end of an energy gap, as 
expected from the Peierls effect. We first illustrate the 
self-consistent potentials which show no surprise. At low 
density, one recognizes the shapes of clearly resolved indi- 
vidual baryons from Ref. [l7|, see Figs. [T] and [51 At high 
density where the baryons overlap significantly, the low- 
est Fourier modes (So, S*!, Pi) dominate, as anticipated 
in our perturbative calculation (Figs. [3] and [3]). Increas- 
ing 7 tends to wash out the oscillations at all densities. 
The energy difference between the crystal and the ho- 
mogeneous phase is shown in Fig. [5] for 3 values of 7. As 
expected, the crystal phase is favored at all densities and 
7 parameters. The horizontal lines at large pf show the 
asymptotic prediction of Eq. ((22l) , whereas the slopes of 
the straight lines near pf — have been obtained from 
the baryon masses [iTl and the mixed phase of the ho- 
mogeneous calculation (see the appendix of [13] )■ This 
provides us with yet another useful test of the compu- 
tations. Fig. [5] shows the ^/-dependence of the energy 
density, now relative to the vacuum, for the same three 
values of 7. The dots are numerical results. The curves 
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FIG. 3: Same as Fig. [T] but at pf = 2.5 where the baryons 
overlap strongly. 
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FIG. 6: Ground state energy density of crystal at T = 
as a function of pf for 7 = 0.2, 0.6, 1.0 from bottom to top. 
Points: numerical HF results, curves; asymptotic predictions 
according to Eqs. p3p and H23p . matched at the point marked 
by a cross. 
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FIG. 4: Same as Fig. [2] but at pf = 2.5. 
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FIG. 5: Difference between energy density of solitonic crystal 
phase and homogeneous phase versus pf for 3 different values 
of 7. The straight line segments drawn show the analytical 
expectations for small and large pf, respectively, see the main 
text. 



have simply been obtained by matching the asymptotic 
expansions Eqs. (|13p . (P5)) . at the point where they coin- 
cide (indicated by the cross). At the scale of the figure, 
the agreement is perfect, reminiscent of similar findings 
in an earlier numerical study of the non-chiral GN model 



III. CONSTRUCTING THE PHASE DIAGRAM 
A. Grand canonical potential 

The phase diagram in the temperature-chemical po- 
tential plane is best analysed via the grand canonical 
potential density ^. The evaluation of ^ in the relativis- 
tic HF approach is well understood and follows earlier 
studies of the non-chiral GN model, the only small com- 
plication being the fact that the spectrum is no longer 
symmetric under E -^ —E. The main building block is 
the familiar single particle contribution to '5, 



* 



1 

']3L 



Z^ 



In 1 



'P{E^,^-li) 



(24) 



For large positive or negative energy eigenvalues, one has 
to use perturbation theory in order to do the renormal- 
ization analytically. The corresponding expression is 



^ port 



iris:'"(i+-'"^'--"') (2^) 



where Eri^p denotes the 2nd order perturbative eigenvalue 
of the Dirac-HF Hamiltonian H. The standard manipu- 



lation 
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1 r^lnfl + e-^(^+---^)) 
P Jp 27r V / 

2 /-""dp 



/Sip 2^ 



Infl + e'^f^-i'"-^)) 



/■^/2 dp 



(26) 



isolates the divergence in the sum over single particle 
energies, which can then be dealt with like at T = 0, 
see Sec. |ll] and Ref. [l3|: adding the double counting 
correction and using the gap equation to eliminate un- 
physical parameters. We then minimize ^ with respect 
to the potentials S', P. The result is the renormalized 
grand canonical potential density, together with the self- 
consistent potential at a given temperature and chemical 
potential. A vacuum subtraction finally normalizes ^ to 
at the point (T = 0, /^ = 0) and removes remaining 
trivial divergences from the Dirac sea. 



B. Ginzburg-Landau theory 

There are regions in (7,/i, r)-space where a full HF 
calculation can be bypassed. This is the case whenever a 
microscopic GL theory can be derived, leading to an ef- 
fective bosonic field theory directly in terms of the scalar 
and pseudoscalar potentials 5*, P with the fermions "inte- 
grated out" . One can identify two such regions requiring 
somewhat different approximations. Close to the tricrit- 
ical point at (7 = 0, ^ = 0,r = e^/n), the potentials 
are both weak and slowly varying. This was exploited 
in Ref. [la], where a GL effective action was obtained 
analytically, using the derivative expansion around the 
free, massless fermion theory. The resulting effective 
action was then minimized by a numerical solution of 
the Euler-Lagrange equation, an inhomogeneous, com- 
plex non-linear Schrodinger equation. In this manner, 
a soliton crystal solution could be identified in a small 
region of (7,/x, T) space, separated by 2nd and 1st or- 
der transitions from a homogeneous massive Fermi gas 
phase. We refer the reader to this paper for more de- 
tails. Another approximation allows one to study the 
phase diagram for 7 <C 1 and ^ <gC 1, but without any 
restriction in temperature. Here, the potentials are still 
slowly varying but develop a large, constant scalar term 
5*0, i.e., a mass. The derivative expansion can still be 
trusted, provided one expands now around the massive 
free Dirac theory. This technique was applied some time 
ago at T = to the baryons in the chiral GN model near 
7 = [l6| . The generalization to finite temperature and 
chemical potential is technically rather involved. In par- 
ticular, it does not lead anymore to analytic expressions 
as in Ref. [lg|, since the thermal integrals with massive 
single particle energies cannot be done in closed form. 
As this technique was used here only for a small part of 
the phase diagram, we refrain from giving all the details 



which have been worked out in Ref. (2l|. The resulting 
effective action is a polynomial in S, P and its deriva- 
tives with (7, /i, T)-dependent coefficients given in terms 
of one-dimensional numerical integrals. It can be mini- 
mized numerically by varying the period and the Fourier 
coefficients of S and P, resulting in the equilibrium value 
of $. In this way, it is possible to extend the calculation 
of the 1st order transition line at small 7 down to zero 
temperature and check that the base point of the criti- 
cal line coincides with the baryon mass. Some examples 
of results for the phase boundary thus obtained will be 
shown below together with the results of the full HF cal- 
culation, see Sec.|niD]and Figs. \U\ fT2l 



C. Perturbative 2nd order phase boundary 

As is well understood by now from similar studies of 
the non-chiral GN model or from the GL approach near 
the tricritical point, the exact location of a contingent 
2nd order phase boundary between crystal and homoge- 
neous phases is a perturbative matter. For this purpose, 
5*0 (i.e., the dynamical fermion mass) has to be treated 
exactly, whereas it is sufficient to keep Si , Pi from the in- 
homogeneous terms and treat them in 2nd order almost 
degenerate perturbation theory (ADPT). As a matter of 
fact, right at the phase boundary this amounts to naive 
2nd order perturbation theory and a principal value pre- 
scription for integrating through the pole when summing 
over single particle states [20] • The Hamiltonian is di- 
vided up according to 



H^Ho + V 



(27) 



where 



Ho 



1 d 
1 ox 



V = -f°2SiCOs{2pfx) ~ i'y^2Pism{2pfx). (28) 

To define the notation, we cast the unperturbed problem 
into the form (ry = ±1 is the sign of the energy) 



HalViP) =11^17], p), E = yjp^ + m^ 

with the free, massive spinors 



ipa; 



Matrix elements of V are then given by 
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with 
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[riE{ip' + m) — i]' E'{[p — m)] 
i [rjri' EE' + (ip — m)(ip' + m)] 

>^{Sp',p+2pf 



^p' .p~2pf )i 



(32) 



leading to the following 2nd order energy shift, 



6E„ 



We insert 



{p^-p))E 



Er),p — rjE + SEri^p 



(33) 



(34) 



into the single particle contribution to the grand canon- 
ical potential density, 



* 



2 f^'^dp 
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(l + e-'^'^-i'"-'^)) 



{l + e-P^'^^^'-^''^) 



(35) 



and linearize in SE^jp. Adding the usual HF double 
counting correction term and invoking the gap equation 
for the fermion mass at finite T, /i in the translationally 
invariant case. 



= ?Ti(7 + Injn) — 7 

dp / 1 1 



(36) 



to simplify the resulting expression, the perturbative cor- 
rection to the grand canonical potential becomes 
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Sl + P^l EJSf+p}P?^fEf-p, 



2TrpfEf 
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Ef+pf 



The energies E, Ef are defined with the mass m = Sq and 
momenta p,Pf, respectively. The principal value inte- 
grals are the only remnant of ADPT at the phase bound- 
ary [20| and have to be evaluated numerically. The phase 
boundary can now be found using the following strategy. 
In 2nd order perturbation theory, according to Eq. (|37p 
we may write the grand canonical potential schematically 
as 

* = *hom + MnSl + 2M12S1P1 + M22P1 (38) 

where all coefficients depend on m and pf. We have to 
vary ^ with respect to the 4 parameters, m,Si,Pi and 



Pf. This yields the 4 equations 



= 



= 
= 

= 



dm 
+2S1P1 

SiMn 
S1M12 

Si 



dM 



dr 



dMi2 

dm 
f P1X12 

f P1M22 
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p^ 



,dM 



22 



dm, 



2dMii 



dpf 



2SiP, 



dMi2 
dpf 



Pi 



dM 



22 



dpf 



(39) 

(40) 
(41) 

(42) 



At the phase boundary, Eq. ([55)1 can be simplified to 
the standard equation for the homogeneous phase since 
^i. Pi vanish. 



dm 



= 0. 



(43) 



Eqs. (|40p . p4|) represent a homogeneous system of equa- 
tions which can be cast into the equivalent form 



detM = MiiM22-Mi2 = 
Sj_ ^ M12 
Pi " Mil' 

Dividing Eq. (|^ by Pf and using Eqs. 
finally obtain the condition 



0, 



ddetM 
dpf 



= 0. 



(44) 
(45) 

(|i5)) . we 
(46) 



In order to determine the phase boundary, we have to 
find the points in the (/x,T) plane where Eqs. (gSl), ((44)) 
and (j46| hold simultaneously. Eq. ([45]) then yields the 
unstable direction. All of this can be done numerically 
to any desired accuracy. Before turning to the results, it 
may be worthwhile to ask whether we can say anything 
about the outcome of the calculation beforehand. Indeed, 
it is easy to determine the asymptotic behaviour of the 
perturbative 2nd order sheet in the limit /i — > 00, for any 
7. Along the lines of a similar analysis in the appendix 
of Ref. [20| we arrive at the approximate expression for 
the grand canonical potential valid at large /i « p/. 



* = 



02 

27r 
,2 



[7 + ln{2pf)] - 2:^ + ^ [7 + ln(4p/)] (47) 



TT 

00 



+ ^(27 - 1 + lny2) _ / dpln(l + e-^Vp^+y^) 

An (JTT J„ 

(using once again variables X — {Si+Pi)/2, y = Si— Pi). 
Sq and X are not affected by finite temperature at all, so 
that Eqs. (flS)) still hold. Minimization with respect to y 
yields either y = (translationally invariant solution) or 
the condition 



7 + In 2/ + 2 



dp 



1 



1 



v^ 



y2 ^(3^/^^^ 



1 



0. (48) 



In order to compute the phase boundary, we expand the 
integral for small y [22|, 



7 + In y — In 



Py 



C + 0(2/2 



0, 



(49) 




FIG. 7: 3d plot of the perturbative, 2nd order phase bound- 
ary in the chiral GN model (upper sheet), compared to the 
corresponding phase boundary in the non-chiral GN model 
(lower sheet). The fat line is the tricritical curve of the latter 
model. The tricritical curve of the chiral GN model cannot 
be determined by this calculation. The dashed curve at /i = 2 
is the asymptotic prediction of Eq. (|50|l . 



where C is the Euler constant. The critical Une where 
the non-trivial solution for y disappears is then given by 
the following asymptotic expression valid at large fi, 



(50) 



Fig. [7] shows the results for the perturbative 2nd or- 
der sheet (a preliminary version of this plot has been 
given before in [2^). This figure actually contains the 
2nd order sheets for both the chiral and the non-chiral 
GN models to highlight the differences between the two 
models. The lower sheet ending at the fat black tricrit- 
ical line belongs to the GN model with discrete chiral 
symmetry. To test our method, we have recalculated 
the curves shown here perturbatively. They agree indeed 
with the results of Ref. 24] where the same critical sur- 
face was deduced from the full, analytical solution of the 
HF problem. The upper sheet in Fig. [7] is the new result 
for the chiral GN model. Here we have supplemented the 
equidistant curves at 7 = 0.1, 0.2, ...2.0 by 2 more curves 
at the small 7 values 0.01 and 0.0001. This is useful to 
illustrate how this 2nd order sheet goes over into the hor- 
izontal critical line T = e /tt in the chiral limit 7^0. 
We also compare the 2nd order sheet with the analyti- 
cal prediction, Eq. ([50)1 . at large fi. As the dashed curve 
shows, the full results are already indistinguishable from 
this formula at /i = 2. At low /i, the curves bend over. 
Our calculation gives us no clue as to where the tricrit- 
ical points are beyond which these curves turn into 1st 
order critical lines. We will get back to this issue in the 
following subsection when we discuss the full HF calcula- 
tion. Notice also that at large 7, the perturbative sheets 
of both variants of the GN model seem to come together 



at the same line (the tricritical line of the discrete chiral 
GN model) , whereas this does not hold anymore at small 

7- 

Finally, we should stress the fact that the sheet in 
Fig. [7] represents the surface where the homogeneous 
phase becomes unstable towards crystallization in a con- 
tinuous transition. If a 1st order transition occurs before 
reaching this sheet from the outside, there will be no 2nd 
order transition and the corresponding part of the 2nd 
order sheet becomes obsolete. As shown below, this is 
indeed what happens at sufficiently low temperatures. 



D. Non-perturbative 1st order phase boundary 
and full phase diagram 

The most tedious task of the present study is to deter- 
mine the 1st order phase boundary. For arbitrary 7 and 
/i, no shortcut like GL theory is known and we have to 
resort to the full, numerical HF calculation. For a given 
point in the (7, /i, T) diagram, we evaluate the renormal- 
ized grand canonical potential density \1/ by minimization 
with respect to the period a and the Fourier components 
Se,Pi of the mean field. The critical line is then con- 
structed as follows. We evaluate ^ along a straight line 
trajectory for fixed 7, T and several equidistant values 
of /i, starting from inside the anticipated crystal phase 
and proceeding towards lower /x values. We then plot 
^ against /x and compare this thermodynamic potential 
with the one of the homogeneous solution. In Fig. [51 we 
illustrate the outcome of such a computation for the case 
7 = 1.0, T = 0.08. The thermodynamically stable phase 
is the one with the lowest value of \1/, hence the point 
of intersection of the 2 curves defines the critical chem- 
ical potential at this temperature. Since we can follow 
the crystal solution beyond this point (before it jumps 
onto the other curve), this is clearly a 1st order transi- 
tion where two different solutions coexist at the phase 
boundary. The difference in slopes at the intersection 
point translates into two different densities, so that a 
mixed phase would appear in a (p, T) phase diagram. 
The critical point can be determined accurately in cases 
like that shown in Fig. [51 By contrast, Fig. [51 illustrates 
an example where the transition is likely to be 2nd order, 
namely at 7 = 1.0, T = 0.12. Here, one does not see a 
crossing of the two curves. Due to the limited numeri- 
cal accuracy, one cannot rule out a very weak first order 
transition, therefore it is difficult to locate the tricritical 
point precisely in this manner. 

The result of such a computation of the 1st order criti- 
cal line at 7 = 1.0 is shown in Fig.[TOl The solid line is the 
perturbative 2nd order line from Sec. piG[) and Fig. [3 
without information on the tricritical point. The squares 
are numerically determined 1st order phase transitions. 
We only show those points for which we could unam- 
biguously identify a 1st order transition. Above T = 0.1, 
there was no visible line crossing anymore. In this way, 
a small gap between the 2nd and 1st order phase bound- 
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FIG. 8: Determination of the 1st order phase boundary at 
T = 0.08,7 = l-O- Points: Grand canonical potential from 
numerical HF calculation vs. /i. The crystal phase can be 
followed down to /i = 0.88. Solid line: Prediction assuming 
homogeneous condensates only. The crossing of the 2 lines 
yields the critical chemical potential for a 1st order transition. 
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FIG. 10: Example of a construction of the phase boundaries 
at 7 = 1.0. Solid line: 2nd order, perturbative critical line 
from Fig. [T] Points: 1st order, non-perturbative critical line 
determined as shown in Fig. [S] The tricritical point has not 
yet been located but must lie on the solid line, above T — 0.1. 
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FIG. 9: Same plot as in Fig. |8] at T = 0.12,7 = 1.0. The 
absence of line crossing is indicative of a continuous, 2nd order 
phase transition. 



aries is left. All we can say is that the tricritical point 
lies on the 2nd order line above the last 1st order point 
shown, i.e., at T > 0.10 in the case at hand. For as much 
as we can tell, the two critical lines are joined tangen- 
tially at the tricritical point. Note that the base point 
of the 1st order line at T = drawn here is the baryon 
mass at 7 = 1.0 taken from Ref. [l3|- The fact that the 
numerical points interpolate nicely between the baryon 
mass at T = and the perturbative phase boundary is a 
healthy sign, suggesting that the accuracy reached here 
is adequate. 

In a lengthy numerical calculation with MAPLE, we 
have determined a number of 1st order critical lines, see 



0.6 ^ 


^ 


■s^<r°*°°^ 


0.5^ 




/\5^ 


0.4: 




f \/ 


T 0.3: 




1 \j 


0.2: 




1 1 


0.1: 
0- 


■c^ 





0.4" 



0.6 



1.2 



1.4 



1.6 



1.2 0.2' 



FIG. 11: Summary of all the results about the phase diagram 
of the chiral GN model obtained in this work. Fat solid curve 
at T = 0: Baryon mass, solid lines at fixed 7: Perturbative 
2nd order sheet, points: numerically determined 1st order 
sheet, computed in steps of A7 = 0.1, AT = 0.01. The 2 
curves at very small 7 are taken from the GL analysis [211 ] 
and belong to 7 = 0.01 and 0.0001, respectively. The fat line 
crossing the 2nd order sheet is the tricritical curve. 



Fig. [TTJ The solid curve at T = is the baryon mass 
from Ref. T7]. The thin lines are the 2nd order critical 
lines from Fig. [T] the points arc numerically determined 
1st order transitions computed on a grid with resolution 
A7 — 0.1, AT = 0.01. Also shown are two additional 
curves at very small 7 (0.01 and 0.0001) obtained pre- 
viously by means of the GL theory f21|. These results 
confirm the picture discussed in connection with Fig. [TU] 
and provide us with a first candidate for the full phase di- 
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FIG. 12: DifFerence in grand canonical potential between 2 
phases near the tricritical point, using GL theory near 7 = 0. 
The rescaled potential difference is plotted vs. a ~ {Tc~tY''^ 
along the 2nd order critical line. 
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FIG. 13: Square of rescaled Fourier amplitudes Si (solid) 
and Pi (dashed) vs. g for the calculation corresponding to 
Fig. [12] The linear behavior shows that Si, Pi - (Tc - Tf^* 
and locates the tricritical point precisely ai a — 1.464. 



agram of the chiral GN model. We find no indication of 
any further phase transitions beyond those which have 
been identified in the earher study near the tricritical 
point [Hi. 

In Fig. [Til we have also plotted a tricritical line where 
the 1st and 2nd order critical sheets are joined together. 
As is clear from the gap between the calculated 1st order 
sheet and this line, some extrapolation had to be used. 
We proceeded as follows. For a fixed value of 7, we move 
along the 2nd order instability line, starting well below 
the expected tricritical point. We then perform the HF 
minimization and follow in particular the evolution of 
the largest Fourier components Si, Pi. At the tricritical 
point, these are expected to vanish with some power law 
~ [Tc — T)". In order to find the relevant critical expo- 
nent a, we went back to the GL approach near 7 = [l8| 
and performed a similar analysis there. This has the ad- 
vantage that we can work with a much higher numerical 



precision in this regime. Let us first recall that to take 
advantage of simple scaling properties near the tricriti- 
cal point, the variables /i, T have been replaced by the 
rescaled variables 



27- 



-1/3 



^J■, 



-1/3 



VTc-T (51) 



with a ^ 6.032, Tc = 0.5669 in Ref. [l8|. We now move 
along the perturbative phase boundary plotted in Fig. 6 
of [la] between cr = 1.4 and 1.6 enclosing the tricritical 
point. Along this trajectory the effective action is min- 
imized with respect to the Fourier components of Se,Pe 
{£ < 4) and the period. The resulting grand canonical 
potential is compared to the homogeneous calculation in 
Fig. [m Fig. [131 then shows clearly that the Fourier com- 
ponents Si, Pi vanish like a^^^ ~ (T^ - Py/^. (Notice 
that the grand canonical potential and the Fourier com- 
ponents in Figs. [T2| and [T3| have been rescaled by the 
factors 27ra/7 and 7^^/^^, respectively, cf. Ref. 18].) As 
a by-product, we have determined in this way a more ac- 
curate value of the tricritical point near 7 = than in 
Ref. [1^, namely at = 1.464, i^t = 3.039. Coming back to 
the full HF calculation, we have located the point where 
Si , Pi vanish along the 2nd order instability curve assum- 
ing the same critical exponent a = 1/4 for all 7. Due to 
numerical limitations, the extrapolation is not as quanti- 
tative as in Fig. [T21 but still fairly straightforward. The 
result is the tricritical curve drawn in Fig. 1111 



IV. SUMMARY AND CONCLUSIONS 

To summarize, we have redrawn the phase diagram of 
Fig. [TT] in a way which shows more clearly the shape 
of the 2 critical sheets, see Fig. [TH Here, we hide the 
"engineering details" of the underlying construction still 
visible in Fig. [11] Whereas the vertical (1st order) lines 
have actually been computed via HF, the horizontal lines 
are composed of straight line segments joining neighbor- 
ing points to guide the eye. The tricritical line separates 
1st and 2nd order sheets. It is not completely smooth 
since this particular curve is the most difficult part of 
the whole calculation, exceptionally sensitive to numeri- 
cal inaccuracies. 

It is interesting to compare this newly determined 
phase diagram of the massive chiral GN model to other 
related phase diagrams. For this purpose, we have taken 
the results for the discrete chiral GN model from Ref. [2^ 
and plotted them at the same scale and under the same 
viewing angle as in Fig. [TJ] cf. Fig. [13 Here the 2 criti- 
cal sheets are both 2nd order and joined in a cusp rather 
than tangentially. The qualitative differences between 
Figs. [Til and [T5I are due to the difference between contin- 
uous and discrete chiral symmetries of the two GN-type 
models, reflecting the corresponding universality classes. 
If one would only admit homogeneous phases as was done 
in the early works on these phase diagrams, the 2 models 
would give identical results. This is illustrated in Fig. [TBI 
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FIG. 14: Phase diagram of the massive chiral GN model. 
The crystal phase with complex order parameter is separated 
from the massive Fermi gas phase by 1st (dark shaded) and 
2nd (light shaded) order critical sheets joined at a tricritical 
line. 




FIG. 16: Common phase diagram for both variants of the 
massive GN model, assuming homogeneous condensates only. 
There is only a single massive Fermi gas phase, but the value 
of the mass changes discontinuously across the dark shaded 
1st order sheet. The tricritical line agrees with the one in 
Fig.[TSl Adapted from Refs. [3 and [Sj. 




FIG. 15: Phase diagram of the massive discrete chiral GN 
model, adapted from Ref. [241 ). The crystal phase with real 
order parameter is separated by two 2nd order critical sheets 
from the massive Fermi gas phase, meeting at the tricritical 
line. 



adapted from Ref. [2j]. Here, the dark shaded sheet is 
1st order, and there is only a single massive Fermi gas 
phase at 7 > 0. 

Let us finally comment on some open questions. As 
pointed out above, to determine the tricritical line of the 
chiral GN model requires some extrapolation of numeri- 
cal results. The minimization becomes difficult close to 
the tricritical line where the effective potential is flat. 
Independent analytical work on the tricritical line and 
the critical exponent a discussed above would therefore 
be useful. One would also like to know the universal- 



ity classes to which the different variants of GN models 
belong. 

Another issue where further work is needed is related to 
the possibility of chiral twist in the massless case, recently 
discovered in Ref. [1J| . From the symmetry point of view, 
the situation in the chiral limit may be characterized as 
follows: The Hamiltonian commutes with the generators 
(P, Q, Q5) of translations and vector/axial vector phase 
transformations of the fermions. A mass term (like in the 
vacuum or any homogeneous phase) breaks Qs, reducing 
chiral symmetry to U(l) vector transformations with the 
appearance of a massless Goldstone boson, the pion, and 
leaves P unbroken. The chiral spiral solution breaks P 
and Qs , but leaves the linear combination P + /iQs un- 
broken. Since one unbroken, continuous symmetry is left, 
one expects only one gapless excitation, a mixture of a 
phonon and a pion. If the twisted kink crystal is realized, 
the symmetry will be further broken down to one dis- 
crete combination of translation and 75 phase rotation, 
cf. Eq. pT|) . Such crystals should feature two different 
gapless excitations, the phonon and the pion. In view of 
these different physics implications, it is important to re- 
consider the phase diagram in the chiral limit once again 
and establish the thermodynamically most stable phases. 

In the massive chiral GN model, chiral symmetry is ex- 
plicitly broken by the bare mass term, and Qs does not 
commute with H anymore. It is therefore unlikely that 
the Qs operator appears in a residual discrete symmetry 
of the condensates, as in the chirally twisted kink crystal. 
The only remaining issue is then the fate of translational 
symmetry and its generator P. So far, we have tacitly 
assumed that translational invariance breaks down to a 
discrete subgroup with concomitant periodic potentials. 
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as is comiTLon in condensed matter systems. If this as- 
sumption would turn out to be wrong, the present cal- 
culation should be regarded as a variational calculation 
rather than the exact solution of the model in the large 
N limit, but such a scenario does not seem very likely to 
us. 
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